{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### From Chapter 1 of Twomey (1977)\n",
    "\n",
    "$$ B(\\lambda, T) = \\alpha_{\\lambda} B(\\bar{\\lambda}, T) + \\beta_{\\lambda} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bt2rad import bt2rad\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = np.arange(230,300,10)\n",
    "B = bt2rad(675., T)            # 675 cm-1 is near the center for the 15-um CO2 band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Radiance')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAZ5ElEQVR4nO3df5RcZZ3n8fen0zQoiDSdhuVXOrQb2QFmRLrNlO6oKA4gM0PQOWgYdicHzGRRdITZcQcOuwIic5jFcZzjgE4MSHR6wKiwoHsUQgblzBwb7I4Rwo9IaGiIxKRtGnSEpdPUd/+4t28qTaW70nTVrer6vM7pU7ee+6O+D0X60/e5Vc9VRGBmZgbQkncBZmZWPxwKZmaWcSiYmVnGoWBmZhmHgpmZZVrzLuC1WLhwYSxevDjvMszMGsrg4OAvI6Kz3LqGDoXFixczMDCQdxlmZg1F0vDe1nn4yMzMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzBrM4PAY19+7lcHhsTk/dkN/T8HMrNkMDo9x3pp+xieKtLW20LeyQE9X+5wd32cKZmYNpH9olPGJIsWAXRNF+odG5/T4DgUzswZS6O6grbWFBYL9WlsodHfM6fE9fGRm1kB6utrpW1mgf2iUQnfHnA4dgUPBzKzh9HS1z3kYTPLwkZmZZaoWCpJukrRT0uaStqslPShpk6S7JR1Zsu4ySVslbZF0erXqMjOzvavmmcLNwBlT2q6LiN+JiJOA7wKfBpB0PLAcOCHd5wZJC6pYm5mZlVG1UIiI+4DnprT9quTpgUCky8uAWyPi5Yh4EtgKLK1WbWZmVl7NLzRLugb4U+AF4D1p81FAf8lm29I2MzOroZpfaI6IyyPiGKAP+HjarHKblttf0ipJA5IGRkZGqlWmmVlTyvPTR/8M/HG6vA04pmTd0cCz5XaKiNUR0RsRvZ2dZW8xamZms1TTUJC0pOTpWcBj6fKdwHJJ+0s6FlgCPFDL2szMrIrXFCTdApwCLJS0DbgCOFPScUARGAYuBIiIhyWtAx4BJoCLIuKVatVmZmblKaLs0H1D6O3tjYGBgbzLMDNrKJIGI6K33Dp/o9nMzDIOBTMzyzgUzMws41AwM7OMQ8HMzDIOBTMzyzgUzMws41AwM7OMQ8HMzDIOBTMzyzgUzMws41AwM7OMQ8HMzDIOBTMzyzgUzMws41AwM7OMQ8HMzDJVCwVJN0naKWlzSdt1kh6T9KCk2yUdUrLuMklbJW2RdHq16jIzs72r5pnCzcAZU9rWAydGxO8APwMuA5B0PLAcOCHd5wZJC6pYm5mZlVG1UIiI+4DnprTdHRET6dN+4Oh0eRlwa0S8HBFPAluBpdWqzczMysvzmsIFwPfS5aOAZ0rWbUvbXkXSKkkDkgZGRkaqXKKZWXPJJRQkXQ5MAH2TTWU2i3L7RsTqiOiNiN7Ozs5qlWhm1pRaa/2CklYAfwicGhGTv/i3AceUbHY08GytazMza3Y1PVOQdAbwV8BZEfFiyao7geWS9pd0LLAEeKCWtZnZ/DY4PMb1925lcHgs71LqWtXOFCTdApwCLJS0DbiC5NNG+wPrJQH0R8SFEfGwpHXAIyTDShdFxCvVqs3Mmsvg8BjnrelnfKJIW2sLfSsL9HS1511WXapaKETEuWWab5xm+2uAa6pVj5k1r/6hUcYnihQDdk0U6R8adSjshb/RbGbzXqG7g7bWFhYI9mttodDdkXdJdavmF5rNzGqtp6udvpUF+odGKXR3+CxhGg4FM2sKPV3tDoMKePjIzMwyDgUzM8s4FMzMLONQMDOzjEPBzMwyDgUzM8s4FMzMLONQMDOzjEPBzMwyDgUzM8s4FMzMLONQMDOzjEPBzMwyVQsFSTdJ2ilpc0nbOZIellSU1Dtl+8skbZW0RdLp1arLzMz2rppnCjcDZ0xp2wx8ELivtFHS8cBy4IR0nxskLahibWZmVkbVQiEi7gOem9L2aERsKbP5MuDWiHg5Ip4EtgJLq1WbmZmVVy/XFI4Cnil5vi1texVJqyQNSBoYGRmpSXFmZs2iXkJBZdqi3IYRsToieiOit7Ozs8plmZk1l3oJhW3AMSXPjwaezakWM7OmVS+hcCewXNL+ko4FlgAP5FyTmVnTaa3WgSXdApwCLJS0DbiC5MLzF4FO4P9K2hQRp0fEw5LWAY8AE8BFEfFKtWozM7PyqhYKEXHuXlbdvpftrwGuqVY9ZmY2s3oZPjIzszrgUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDCzvRocHuP6e7cyODyWdylWI1Wb5sLMGtvg8BjnrelnfKJIW2sLfSsL9HS1512WVZnPFMysrP6hUcYnihQDdk0U6R8azbskqwGHgpmVVejuoK21hQWC/VpbKHR35F2S1YCHj8ysrJ6udvpWFugfGqXQ3eGhoyZRUShIOhz4a+DIiHi/pOOBt0fEjVWtzsxy1dPV7jBoMpUOH90M3AUcmT7/GXBxNQoyM7P8VBoKCyNiHVAEiIgJYNo7o0m6SdJOSZtL2g6VtF7S4+lje8m6yyRtlbRF0umz6IuZmb1GlYbCbyR1AAEgqQC8MMM+NwNnTGm7FNgQEUuADelz0uGo5cAJ6T43SFpQYW1mZjZHKg2FvwDuBN4k6d+ArwGfmG6HiLiP5J7MpZYBa9PltcDZJe23RsTLEfEksBVYWmFtZmY2Ryq60BwRGyW9GzgOELAlInbN4vUOj4jt6TG3SzosbT8K6C/Zblva9iqSVgGrABYtWjSLEszMbG8qOlOQdBFwUEQ8HBGbgYMkfWwO61CZtii3YUSsjojeiOjt7OycwxLMzKzS4aM/i4jnJ59ExBjwZ7N4vR2SjgBIH3em7duAY0q2Oxp4dhbHNzOz16DSUGiRlP01n14EbpvF690JrEiXVwB3lLQvl7S/pGOBJcADszi+mZm9BpV+o/kuYJ2kL5MM61wIfH+6HSTdApwCLJS0DbgCuDY9zkeAp4FzACLiYUnrgEeACeCiiJj2I69mZjb3FFF26H7PjaQW4L8Bp5KM/98NrMn7F3dvb28MDAzkWYKZWcORNBgRveXWVfrpoyLwpfTHzMzmqUrnPvrPwJVAV7qPgIiI7uqVZmZmtVbpNYUbgUuAQWaY3sLMzBpXpaHwQkR8r6qVmJlZ7ioNhXslXQfcBrw82RgRG6tSlZmZ5aLSUPjd9LH0anUA753bcszMLE+VfvroPdUuxMzM8lfx7Tgl/QHJ1NYHTLZFxGeqUZSZmeWj0gnxvgx8mGS6bJF8E7mrinWZmVkOKp376B0R8afAWERcBbydPSewMzOzeaDSUHgpfXxR0pHALuDY6pRkZmZ5qfSawnclHQJcB2wk+eTRmqpVZWZmuaj000dXp4vflvRd4ICImOkezWZm1mCmDQVJ742If5H0wTLriIjbqleamZnV2kxnCu8G/gX4ozLrguQbzmZmNk9MGwoRcUX6eH5tyjEzszzNNHz0F9Otj4jPz+ZFJX2S5B7PAr4SEV+QdCjwDWAx8BTwofRe0GZmViMzfST1DelPL/BR4Kj050Lg+Nm8oKQTSQJhKfAW4A8lLQEuBTZExBJgQ/rczMxqaKbho6sAJN0NnBwRv06fXwl8c5av+VtAf0S8mB7rh8AHgGUk93QGWAv8APirWb6GWa4Gh8foHxql0N1BT1d73uWYVazS7yksAsZLno+TDPPMxmbgGkkdJF+KOxMYAA6PiO0AEbFd0mHldpa0ClgFsGjRolmWYFY9g8NjnLemn/GJIm2tLfStLDgYrGFUGgpfBx6QdDvJp44+AHxtNi8YEY9K+htgPfDvwE+BiX3YfzWwGqC3tzdmU4NZNfUPjTI+UaQYsGuiSP/QqEPBGkZF01xExDXABcAY8DxwfkT89WxfNCJujIiTI+JdwHPA48AOSUcApI87Z3t8szwVujtoa21hgWC/1hYK3R15l2RWsYqnzo6IQUnPkE6dLWlRRDw9mxeVdFhE7JS0CPggyQR7xwIrgGvTxztmc2yzvPV0tdO3suBrCtaQKgoFSWcBfwscSfIX/CLgMZL7K8zGt9NrCruAiyJiTNK1wDpJHwGeJpme26wh9XS1OwysIVV6pnA1UADuiYi3SnoPcO5sXzQi3lmmbRQ4dbbHNDOz167SqbN3pb+0WyS1RMS9wElVrMvMzHJQ6ZnC85IOAu4D+iTtZB8+MWRmZo2h0jOFZcCLwCXA94EnKD9JnpmZNbBK76fwm3SxCKyVtABYDvRVqzAzM6u9ac8UJB0s6TJJ/yDpNCU+DgwBH6pNiWZmVisznSl8neQLaz8CVgKfAtqAZRGxqcq1mZlZjc0UCt0R8dsAktYAvwQWTU6MZ2Zm88tMF5p3TS5ExCvAkw4EM7P5a6YzhbdI+lW6LOB16XMBEREHV7U6MzOrqZnup7CgVoWYmVn+Kv2egpmZNQGHgpmZZRwKZmaWcSiYmVnGoWBmZhmHgpmZZXIJBUmXSHpY0mZJt0g6QNKhktZLejx99G2rzMxqrOahIOko4M+B3og4EZiccfVSYENELAE2pM/NzKyG8ho+aiX5dnQr8HrgWZJ7NqxN168Fzs6pNjOzplXzUIiInwOfA54GtgMvRMTdwOERsT3dZjtwWLn9Ja2SNCBpYGRkpFZlm5k1hTyGj9pJzgqOBY4EDpT0XyrdPyJWR0RvRPR2dnZWq0wzs6aUx/DR+0hmWx2JiF3AbcA7gB2SjgBIH3fmUJuZWVPLIxSeBgqSXi9JwKnAo8CdwIp0mxXAHTnUZmbW1Cq6R/Ncioj7JX0L2AhMAD8BVgMHAeskfYQkOM6pdW1mZs2u5qEAEBFXAFdMaX6Z5KzBzMxy4m80W90YHB7j+nu3Mjg8lncpZk0rlzMFs6kGh8c4b00/4xNF2lpb6FtZoKfLX2o3qzWfKVhd6B8aZXyiSDFg10SR/qHRvEsya0oOBasLhe4O2lpbWCDYr7WFQndH3iWZNSUPH1ld6Olqp29lgf6hUQrdHR46MsuJQ8HqRk9Xu8PALGcePjIzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4xDwczMMjUPBUnHSdpU8vMrSRdLOlTSekmPp4+e78DMrMZqHgoRsSUiToqIk4Ae4EXgduBSYENELAE2pM/NzKyG8h4+OhV4IiKGgWXA2rR9LXB2blWZmTWpvENhOXBLunx4RGwHSB8PK7eDpFWSBiQNjIyM1KhMM7PmkFsoSGoDzgK+uS/7RcTqiOiNiN7Ozs7qFGdm1qTyPFN4P7AxInakz3dIOgIgfdyZW2VmZk0qz1A4l91DRwB3AivS5RXAHTWvyMysyeUSCpJeD/w+cFtJ87XA70t6PF13bR61mZk1s1xuxxkRLwIdU9pGST6NZGZmOcn700dmZlZHHApmZpZxKJiZWcahYGZmGYdCgxscHuP6e7cyODyWdylmNg/k8ukjmxuDw2Oct6af8Ykiba0t9K0s0NPlyWXNbPZ8ptDA+odGGZ8oUgzYNVGkf2g075LMrME5FBpYobuDttYWFgj2a22h0N0x805mZtPw8FED6+lqp29lgf6hUQrdHR46MrPXzKHQ4Hq62h0GZjZnPHxkZmYZh4KZmWUcCmZmlnEomJlZxqFgZmYZh4KZmWXyuvPaIZK+JekxSY9KerukQyWtl/R4+ujPWZqZ1VheZwp/D3w/Iv4T8BbgUeBSYENELAE2pM/NzKyGah4Kkg4G3gXcCBAR4xHxPLAMWJtuthY4u9a1mZk1uzzOFLqBEeCrkn4iaY2kA4HDI2I7QPp4WLmdJa2SNCBpYGRkpHZVm5k1gTxCoRU4GfhSRLwV+A37MFQUEasjojciejs7O6tVo5lZU8ojFLYB2yLi/vT5t0hCYoekIwDSx5051GZm1tRqHgoR8QvgGUnHpU2nAo8AdwIr0rYVwB21rs3MrNnlNUvqJ4A+SW3AEHA+SUCtk/QR4GngnJxqMzNrWrmEQkRsAnrLrDq11rWYmdlu/kazmZllHApmZpZxKJiZWcahYGZmGYeCmZllHApmZpZp2lAYHB7j+nu3Mjg8lncpZmZ1I68vr+VqcHiM89b0Mz5RpK21hb6VBXq6fPsGM7OmPFPoHxplfKJIMWDXRJH+odG8SzIzqwtNGQqF7g7aWltYINivtYVCd0feJZmZ1YWmHD7q6Wqnb2WB/qFRCt0dHjoyM0s1ZShAEgwOAzOzPTXl8JGZmZXnUDAzs4xDwczMMg4FMzPLOBTMzCzjUDAzs4wiIu8aZk3SCDD8Gg6xEPjlHJWTp/nSD3Bf6tF86Qe4L5O6IqKz3IqGDoXXStJARJS7V3RDmS/9APelHs2XfoD7UgkPH5mZWcahYGZmmWYPhdV5FzBH5ks/wH2pR/OlH+C+zKiprymYmdmemv1MwczMSjgUzMwsM29DQdIxku6V9KikhyV9Mm2/WtKDkjZJulvSkSX7XCZpq6Qtkk7Pr/rd9taPkvV/KSkkLSxpq7t+wLTvyZWSfp6+J5sknVmyT0P1JV33ibTehyX975L2huqLpG+UvCdPSdpUsk/d9WWafpwkqT/tx4CkpSX71F0/YNq+vEXSjyQ9JOk7kg4u2Wdu+hIR8/IHOAI4OV1+A/Az4Hjg4JJt/hz4crp8PPBTYH/gWOAJYEG99iN9fgxwF8kX+BbWcz9meE+uBP6yzPaN2Jf3APcA+6frDmvUvkzZ5m+BT9dzX6Z5T+4G3p+2nwn8oJ77MUNffgy8O22/ALh6rvsyb88UImJ7RGxMl38NPAocFRG/KtnsQGDySvsy4NaIeDkingS2AkvJ2d76ka7+O+B/sLsPUKf9gBn7Uk4j9uWjwLUR8XK6bme6SyP2BQBJAj4E3JI21WVfpulHAJN/Ub8ReDZdrst+wLR9OQ64L91sPfDH6fKc9WXehkIpSYuBtwL3p8+vkfQMcB7w6XSzo4BnSnbbxvS/sGqutB+SzgJ+HhE/nbJZ3fcDXv2eAB9Ph/VukjR5S7xG7MubgXdKul/SDyW9Ld2sEfsy6Z3Ajoh4PH1e932Z0o+LgevSf/OfAy5LN6v7fsCr+rIZOCtddQ7JaAHMYV/mfShIOgj4NnDx5FlCRFweEccAfcDHJzcts3vdfF63tB/ABHA5uwNtj03LtNVNP6Dse/Il4E3AScB2kqEKaMy+tALtQAH4FLAu/Uu7Efsy6Vx2nyVAnfelTD8+ClyS/pu/BLhxctMyu9dNP6BsXy4ALpI0SDKsND65aZndZ9WXeR0KkvYj+Q/aFxG3ldnkn9l9+rWN3akLcDS7TzNzVaYfbyIZN/yppKdIat0o6T9Qx/2A8u9JROyIiFciogh8hd2nvQ3XF5Kab4vEA0CRZOKyRuwLklqBDwLfKNm8bvuyl36sACaXv0kD//8VEY9FxGkR0UMS1E+km89dX/K+oFKtH5Lk/BrwhSntS0qWPwF8K10+gT0v1AxRBxed9taPKds8xe4LzXXZjxnekyNKli8hGRtt1L5cCHwmXX4zySm9GrEv6bozgB9OaavLvkzznjwKnJIunwoM1nM/ZujL5AcXWtL1F8x1X3LvfBX/o/4eyenTg8Cm9OdMkuTdnLZ/h+Ti8+Q+l5Mk7xbSTyvk/bO3fkzZJguFeu3HDO/J14GH0vY7p4REo/WlDfin9P+xjcB7G7Uv6bqbgQvL7FN3fZnmPfk9YDD9pXk/0FPP/ZihL58k+STSz4BrSWelmMu+eJoLMzPLzOtrCmZmtm8cCmZmlnEomJlZxqFgZmYZh4KZmWUcCtaQJHWUzOD5iymzrLblXV85ki5Iv2BYreMfKOkHklok/ccps5peKOnHkt4o6QuS3lWtOqyxteZdgNlsRMQoybQYSLoS+PeI+FyuRSW1LIiIV/ay+gKS7y78Yh+O1xoRExVuvhL4ZkQUk5k1smOcT/KluvdGxAuSvgj8A7snVjPL+EzB5h1JKyQ9kJ413JD+5dwq6XlJ10naKOkuSb+bTlo3pPQeDpJWSro9Xb9F0v+s8LiflfQAsFTSVelf5ZslfVmJD5OE2OQ9CtokbZN0SHrsgqR70uXPSvpHSeuBr6av8fn0tR+UtHIvXT8PuGPKf4s/Af47cFpEPAcQEU8AR0jqnNP/8DYvOBRsXpF0IvAB4B0RcRLJ2fDydPUbgbsj4mSSicSuJJn24BzgMyWHWZruczLwJ0pu0jLTcTdGxNKI+BHw9xHxNuC303VnRMQ3SL6V+uGIOCkixpneW4E/ioj/CqwCdkbEUuBtJBOiLZrS7wOAoyNiW0lzN/B5kkDYyZ5+ArxjhhqsCXn4yOab95H84hxIh1Bex+4phV+KiPXp8kPACxExIekhYHHJMe6KiDEASf+HZMqB1mmOOw7cXrL/qZI+BRxAMhneIPC9fezHHRHx/9Ll04DfklQaQkuAp0u2Pwx4bsoxdgC/Jpn08YtT1u0EjsRsCoeCzTcCboqI/7VHYzLbZ+lf50Xg5ZLl0n8LU+d+iRmO+1Kk88VIej3JeP3JEfFzSZ8lCYdyJth9tj51m99M6dPHImLDXo4D8NJejvF+4F8l7UzPViYdkO5jtgcPH9l8cw/wIaX3rE4/pbRohn2mOk3SIekv+GXAv+3DcV9HEjK/lPQGdk/NDslf7W8oef4U0JMul2431V3Ax9IAQtJxkl5XukFEjAAHTP3kVUTsIJnp9DpJ7ytZ9WaSSfvM9uAzBZtXIuIhSVcB90hqAXaRfPJmX+aW/1eSe228Cfh6RGwCqOS4ETEqaS3JL9xh9ryD2VeBNZJeIrlucSXwFUm/AB6Ypp5/BBYBm9Khq50kYTXVBpLrBD+YUtMTks4GviNpGcnQ2WKS6wpme/AsqWYl0k/2nBgRF+ddy75ScuvPj0XE+TNsdw5wfERcVZvKrJF4+MhsnoiIH5NcP5jp37WAv6tBSdaAfKZgZmYZnymYmVnGoWBmZhmHgpmZZRwKZmaWcSiYmVnm/wOd26j98bo9bgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(T,B,'.')\n",
    "xlabel('Temperature (K)')\n",
    "ylabel('Radiance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Radiance')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deXxVxf3/8dcnBlxwQ4gIyhbEBVyQRIxb3fdWbC0WSy11KVXRKlWrlLq2+OPrVq2KSnGhLVXRqmitFcSFao2UICCLFIwiCAUaUFAUCPn8/piTcI0hCSH3nru8n49HHvfcOefcfIZAPszMmRlzd0RERADy4g5ARETSh5KCiIjUUFIQEZEaSgoiIlJDSUFERGrkxx3A1mjbtq136dIl7jBERDJKWVnZ/9y9oK5zGZ0UunTpwtSpU+MOQ0Qko5jZws2dU/eRiIjUUFIQEZEaSgoiIlJDSUFERGooKYiISA0lBRERqaGkICKSYcoWruL+1xZQtnBVs392Rs9TEBHJNWULVzFgdCnrK6tomZ/H2ItKKOrcutk+Xy0FEZEMUlpewfrKKqocNlRWUVpe0ayfr6QgIpJBSgrb0DI/j20MWuTnUVLYplk/X91HIiIZpKhza8ZeVEJpeQUlhW2atesIlBRERDJOUefWzZ4Mqqn7SEREaiQtKZjZI2a23MxmJZT9xsxmmtl0M5tgZh0Szg01swVmNs/MTklWXCIiGW3NGhg/Pmkfn8yWwmPAqbXKbnf3g9y9F/A34AYAM+sB9Ad6RveMNLNtkhibiEjmWbMGDjgAvv99+OSTpHyLpCUFd58MrKxVtjrhbSvAo+O+wBPuvs7dPwQWAH2SFZuISEb5/PPwutNOcOmlMHky7LlnUr5VyscUzGy4mS0CBhC1FIA9gUUJly2Oyuq6f5CZTTWzqStWrEhusCIicaqqgpEjoWNHmDIllF17LRx+eNK+ZcqTgrsPc/eOwFjgsqjY6rp0M/ePcvdidy8uKKhzNzkRkcz33ntw5JEweDAUFUGb5p2PsDlxPn30F+Ds6Hgx0DHh3F7AkpRHJCKSDn7zG+jdG+bPhz/+ESZOhG7dUvKtU5oUzKx7wtszgfej4+eB/ma2rZl1BboDU1IZm4hI2mjRAgYMgPffh/POA6urMyU5kjZ5zcweB44F2prZYuBG4HQz2xeoAhYCFwO4+2wzGwfMASqBwe6+MVmxiYiklRUr4KqroG9fOPvsMG6QwkSQKGlJwd3PraP44XquHw4MT1Y8IiJpxz10D111FaxeDYccEspjSgigGc0iIvGYPx9OPBF+8hPYd194910YMiTuqLT2kYhILKZOhbIyePBB+OlPIS89/o+upCAikir/+heUl8OPfgT9+8NJJ0HbtnFH9TXpkZpERLLZp5/CJZfAUUfBrbdCZWUYN0izhABKCiIiyeMOTz8NPXrAqFFw5ZVhZnJ++nbSpG9kIiKZbs4c6NcvPFX0wgthZnKaU0tBRKQ5bdwIr78ejnv2hAkTQusgAxICKCmIiDSfd9+FkhI4/vjQSoAwmJzG3UW1KSmIiGytL76Aq6+GQw+FRYvg8cdh//3jjqpJMid9iYiko8rKkAzmzoVBg2DECGidnP2TU0FJQUSkKVauDL/88/Phmmuge/fwyGmGU/eRiMiWqKqCP/whLGX9zDOh7PzzsyIhgJKCiEjjzZ0LxxwTuokOPjjsl5xllBRERBrj7rtDIpg9Gx5+GF57LSxkl2U0piAiUh/3sCTFHnuEiWi/+x3svnvcUSWNWgoiInVZuRIuvDAkAQgL2I0dm9UJAZKYFMzsETNbbmazEspuN7P3zWymmT1rZrsmnBtqZgvMbJ6ZnZKsuERE6uUOf/kL7LcfjBkDn38ed0QplcyWwmPAqbXKJgIHuPtBwH+AoQBm1gPoD/SM7hlpZtskMTYRkW/68EM49dSwP3LXrmG/gxtuiDuqlEpaUnD3ycDKWmUT3L0yelsK7BUd9wWecPd17v4hsADok6zYRETqtGQJlJbCvfeGvQ8OPjjuiFIuzoHmC4Ano+M9CUmi2uKo7BvMbBAwCKBTp07JjE9EcsGUKfDmm/CLX8CRR8LHH8Muu8QdVWxiGWg2s2FAJTC2uqiOy7yue919lLsXu3txQUFBskIUkWy3Zg38/OdhAbu77w7vIacTAsSQFMxsIPBtYIC7V//iXwx0TLhsL2BJqmMTkRzx3HNhwbr77oPBg2HWLNhpp7ijSgsp7T4ys1OBa4Fj3H1twqnngb+Y2V1AB6A7MCWVsYlIditbuIrS8gqO2rGSg889N6xV9PTToaUgNZKWFMzsceBYoK2ZLQZuJDxttC0w0cwASt39YnefbWbjgDmEbqXB7r4xWbGJSG4pK/8fo4fey8tdD+XeFtvw3F/Gs9+3j4MWLeIOLe0kLSm4+7l1FD9cz/XDgeHJikdEctTMmezZ/8c8MHcG5/a/lSldDmLSrt3ZTwmhTprRLCLZae1auO46KCqizbLFXNP3aqZ0PpAW+XmUFLaJO7q0pbWPRCT7uIdtMP/1Lzj/fFrcfjv9P8+jS3kFJYVtKOqcuZvgJJuSgohkjxUrNm1886tfwQ47wHHHAVDUBiWDRlD3kYhkPnd49NGwXtE994SyM86oSQjSeEoKIpLZ/vMfOOEEuOCCMPfgtNPijiijKSmISOZ6+GE46CCYNg0efBAmT4YePeKOKqNpTEFEMk/1xjf77Qd9+4ZlKtq3jzuqrKCkICKZ49NPw2OmLVvC738fFrA78si4o8oq6j4SkfTnDk89FcYM/vAH2G67UCbNTi0FEUlvixfDxRfDiy9CUVF47d077qiylloKIpLe1q0Lex7cdVfYAEcJIanUUhCR9DNtGjz5JIwYAd26wUcfhYloknRqKYhI+vj8c7jqKjj0UPjjH2Hp0lCuhJAySgoikh5efBF69gzdRD/9KcydCx06xB1VzlH3kYjE74sv4MILoU0b+Oc/4aij4o4oZyWtpWBmj5jZcjOblVDWz8xmm1mVmRXXun6omS0ws3lmdkqy4hKRNFFVBU88AZWV0KoVvPIKvPuuEkLMktl99Bhwaq2yWcD3gMmJhWbWA+gP9IzuGWlm2yQxNhGJ05w58K1vwbnnwrhxoeyAA8KkNIlV0pKCu08GVtYqm+vu8+q4vC/whLuvc/cPgQVAn2TFJiIx+eoruP566NUrjBk8+mhIDJI20mVMYU+gNOH94qjsG8xsEDAIoFOnTsmPTESaz7nnwnPPwY9+FAaUCwrijkhqSZenj6yOsjrnsLv7KHcvdvfiAv2FEkl/FRWwenU4HjoUJkyAP/1JCSFNpUtSWAx0THi/F7AkplhEpDm4w9ixYb2iYcNCWZ8+YZtMSVvpkhSeB/qb2bZm1hXoDkyJOSYRaarycjj11NBNVFgY5h1IRkjamIKZPQ4cC7Q1s8XAjYSB53uBAuBFM5vu7qe4+2wzGwfMASqBwe6+MVmxiUgSPf00/PjHYZ/k++4Li9lto4cJM0XSkoK7b+6Rgmc3c/1wYHiy4hGRJNu4Mfzy7907bHxzxx2wZ53Pi0gaS5fuIxHJVKtXw2WXwXe/G8YRCgvh8ceVEDKUkoKINN1zz4U9kUeOhC5dwuxkyWhKCiKy5ZYtCy2D7343rFf09tthe8wWLeKOTLaSkoKIbLmWLWHGDPi//4OpU+Gww+KOSJpJusxoFpF0N2MG3H03jBoFrVvD++9rraIspJaCiNRv7Vq49tpN+yPPnx/KlRCykpKCiGzehAlh9dLbboOBA0ProEePuKOSJFL3kYjUbePGsDVmixbw+utwzDFxRyQpoKQgIpu4w5//DGeeCbvswnsP/Im3vmjBoV06UBR3bJIS6j4SkWDePDj++LBExejRlC1cRb+Xl3Lb6x8xYHQpZQtXxR2hpICSgkiuW7cObrkFDjoIpk8PTxcNGUJpeQXrK6uocthQWUVpeUXckUoKqPtIJNddcQU89BD84AfhkdM99gCgpLANLfPz2FBZRYv8PEoK28QcqKSCude5l01GKC4u9qlTp8YdhkjmWbUqbI3Zvn1Y5nrePDjttG9cVrZwFaXlFZQUtqGoc+sYApVkMLMydy+u65xaCiK5xB3GjQutg8MOg/HjwwJ2hYV1Xl7UubWSQY7RmIJIrli4EL79bejfH/baC266Ke6IJA0lLSmY2SNmttzMZiWU7WZmE81sfvTaOuHcUDNbYGbzzOyUZMUlkpNeeSVMOnvjjTBu8M47cMghcUclaahRScHM2pnZw2b2UvS+h5ld2MBtjwGn1iq7Dpjk7t2BSdF7zKwH0B/oGd0z0sy0VZPI1lq/PrwWFUG/fjBnTug60k5oshmNbSk8BrwMdIje/we4sr4b3H0yYfvNRH2BMdHxGOCshPIn3H2du38ILAD6NDI2Eant889hyBA4/HDYsCEsYPfYY9CpU9yRSZprbFJo6+7jgCoAd68EmrKHcjt3Xxp9xlJg96h8T2BRwnWLo7JvMLNBZjbVzKauWLGiCSGIZLm//Q169gzdRH36bGotiDRCY5PCF2bWBnAAMysBPmvGOKyOsjqflXX3Ue5e7O7FBQUFzRiCSIb79FM45xz4zndgxx3hzTfhgQegVau4I5MM0thHUn8BPA90M7O3gALg+034fsvMrL27LzWz9sDyqHwx0DHhur2AJU34fJHctcMO8OGH8NvfwjXXaGlraZJGtRTcfRpwDHAE8DOgp7vPbML3ex4YGB0PBMYnlPc3s23NrCvQHZjShM8XyS2zZ8PZZ8Nnn4UkUFoKw4YpIUiTNfbpo8HAju4+291nATua2aUN3PM48Dawr5ktjp5WGgGcZGbzgZOi97j7bGAcMAf4BzDY3ZsyZiGSG776Cn796/BY6RtvhKeKQE8VyVZr1DIXZjbd3XvVKnvX3WN90FnLXEhOevVV+NnPYMGCsKLpnXdC27ZxRyUZpDmWucgzM/Mog0RzCNQ+FYnDbbeF5SomToQTT4w7GskyjU0KLwPjzOxBwlNBFxO6eUQk2dxh7Fg4+mjo3BnGjIGdd4btt487MslCjX0k9VrgVeASYDBhNvIvkxWUiEQ++ABOOQXOOw/uvz+UtWunhCBJ06iWgrtXAQ9EXyKSbBs2wB13hM1vWrYMCeHii+OOSnJAo5KCmR0J3AR0ju4xwN297vV2RWTr3HprWMX07LPh97+HDh0avEWkOTR2TOFhYAhQRtOWtxCRhqxeDcuXw957w89/Dr17h9nJIinU2DGFz9z9JXdf7u4V1V9JjUwklzz7LOy/f1imwj0sYKeEIDFobFJ4zcxuN7PDzax39VdSIxPJBYsXw1lnwfe+BwUF8OCDYHUtBSaSGo3tPjosek2c7ODA8c0bjkgOmToVjjsONm4Mcw+uvBJatIg7KslxjX366LhkByKSM9auDYvXHXwwDBwIV10FXbvGHZUI0PiWAmZ2BmFntO2qy9z9lmQEJZKV1q4NTxQ9+STMnAm77AL33Rd3VCJf09gF8R4EfgBcTngctR/h8VQRaYyXX4YDDoDbb4eTTgqDySJpqLEDzUe4+4+BVe5+M3A4X9//QETq8uWXMGAAnHpqmIT2xhswejTsumvckYnUqbFJ4cvoda2ZdQA2AOoEFWnIdtvBmjVw440wYwZ861txRyRSr8Ymhb+Z2a7A7cA04CPgiWQFJZLR5s2D00+Hjz4Kj5eOHx/GErbdNu7IRBrU2J3XfuPun7r7XwljCfu5+/VN/aZmdoWZzTKz2WZ2ZVS2m5lNNLP50Wvrpn6+SCzWrYObb4aDDoK33w7JATTvQDJKvU8fmdnx7v6qmX2vjnO4+zNb+g3N7ADgp0AfYD3wDzN7MSqb5O4jzOw64DrC6qwi6W/y5LDxzfvvQ//+8LvfwR57xB2VyBZr6JHUYwhLZtc1396BLU4KwP5AqbuvBTCzN4DvAn2BY6NrxgCvo6QgmWLs2LBF5ksvhUFlkQzVqO04m/Ubmu0PjCc8wfQlYW+GqcB57r5rwnWr3L3eLiRtxymxcQ/zDbp1g0MPDYvZbbMNtGoFQNnCVZSWV1BS2IaizuoJlfTS5O04zewX9Z1397u2NBh3n2tm/wdMBD4HZgCVjb3fzAYBgwA6deq0pd9eZOt99BFcemloFfzkJ/Doo2EntEjZwlUMGF3K+soqWubnMfaiEiUGyRgNDTTvFH0VE3Zd2zP6uhjo0dRv6u4Pu3tvd/8WsBKYDywzs/YA0evyzdw7yt2L3b24oKCgqSGIbLnKyrDxTc+e8M9/wj33hDkHtZSWV7C+sooqhw2VVZSWa0FhyRz1thSiiWqY2QSgt7uvid7fBDzV1G9qZru7+3Iz6wR8j9CV1BUYCIyIXsc39fNFkuKRR+Caa+DMM8PyFB3rnr9ZUtiGlvl5bKisokV+HiWFbVIcqEjTNXbto06EJ4WqrQe6bMX3/auZtSFMghvs7qvMbAQwzswuBD4mLKUhEq81a2DBAjjkEDj/fNhrLzjttHofMy3q3JqxF5VoTEEyUmOTwp+AKWb2LOGpo+8Cf2zqN3X3o+soqwBOaOpnijS7F16AwYOhqgo++CBMPjv99EbdWtS5tZKBZKTGTl4bDlwArAI+Bc5391uTGZhIbJYsge9/P3QT7bILPPWUZiNLzmj00tnuXmZmi4iWzjazTu7+cdIiE4nDBx+EvZHXrYNbbw17HbRsGXdUIinTqKRgZmcCdwIdCE8FdQLeJ+yvIJL5PvsstAoKC+GKK+DHP4a99447KpGUa+yCeL8BSoD/uHtX4ETgraRFJZIqX34Jw4ZBly6bFrC75RYlBMlZjU0KG6KB4Dwzy3P314BeSYxLJPkmTQqL1916axg/2HHHuCMSiV1jxxQ+NbMdgcnAWDNbzhbMQhZJK1VVcOGF8NhjoUXwyitwgh58E4HGtxT6AmuBIcA/gA+oe5E8kfSXlxdaBb/6VdgrWQlBpEZjH0n9wt2r3L3S3ccA9wNaClIyx/z5cPLJMGVKeP/738Pw4bD99vHGJZJm6k0KZrazmQ01s/vM7GQLLgPKgXNSE6LIVli/PowZHHggvPMOLFoUyrXxjUidGhpT+BNhwtrbwEXANUBLoK+7T09ybCJb5+23YdAgmDUrTEa75x7o0CHuqETSWkNJodDdDwQws9HA/4BO1QvjiaS1118P8w+efx6+oyEwkcZoaExhQ/WBu28EPlRCkLTlDn/9K7z4Ynh/9dUwe7YSgsgWaCgpHGxmq6OvNcBB1cdmtjoVAYo0yqJF0Ldv6CYaOTKUtWgBO+0Ub1wiGaah/RS2SVUgIk2ycWPY22DYsNBSuOOOsEyFiDRJoxfEE0lLL70EV14Z9jgYOTIsVyEiTdbYyWsi6eOLL+CNN8LxGWeE5SpefFEJQaQZxJIUzGyImc02s1lm9riZbWdmu5nZRDObH71qhxL5ppdeggMOCMlg5cow3+D44zXvQKSZpDwpmNmewM+BYnc/ANgG6A9cB0xy9+7ApOi9SLBsGZx7btj5bLvtQnLYbbe4oxLJOnF1H+UD25tZPrADsISwvtKY6PwY4KyYYpN0s3Il9OgBzzwDN98M06fD0d/Y0VVEmkHKB5rd/RMzuwP4GPgSmODuE8ysnbsvja5Zama713W/mQ0CBgF06tQpVWFLHP73P2jbNrQIbroprF20775xRyWS1eLoPmpNaBV0Jezk1srMftTY+919lLsXu3txQUFBssKUOH31Fdx4I3TsuGkBu8svV0IQSYE4Hkk9kTAzegWAmT0DHAEsM7P2USuhPWHbT8k1b7wBP/sZzJsHP/whdO4cd0QiOSWOMYWPgRIz28HMDDgBmAs8DwyMrhkIjI8hNonT5ZfDsceGlU3/8Q8YOxbatYs7KpGcEseYwjtm9jQwjbB727vAKGBHYJyZXUhIHP1SHZvEwD28moVWwS9/CTfcAK1axRuXSI4yr/5HmYGKi4t96tSpcYchTfXhh3DJJXDRRWHNIhFJCTMrc/fius5pRrOk3oYNcPvt0LMnvPVWmKEMlC1cxf2vLaBs4aqYAxTJXVr7SFJr6tTQMpgxI6xqeu+90LEjZQtXMWB0Kesrq2iZn8fYi0oo6qxJ7SKpppaCpNb8+bBiRZiI9txz4bFToLS8gvWVVVQ5bKisorS8IuZARXKTWgqSfOPHh4loF14I/fuHTW923PFrl5QUtqFlfh4bKqtokZ9HSWGbmIIVyW1KCpI8n3wSHjN99lkoKYHzz4e8vG8kBICizq0Ze1EJpeUVlBS2UdeRSEyUFKT5bdwIDz4IQ4eGQeX/9//gqqtCQqhHUefWSgYiMVNSkOY3YwZcdhmceGJIDt26xR2RiDSSBpqleXz5JbzwQjju3RtKS2HCBCUEkQyjpCBb75VX4MAD4ayzoLw8lB12mDa+EclASgrSdCtWwHnnwUknhfGCiROhsDDuqERkK2hMQZpm3brQTbRsGfz61zBsWNgRTUQympKCbJnFi2HPPWHbbeG22+Cgg8JyFSKSFdR9JI2zfj0MHw577x1mI0PYM1kJQSSrqKUgDfvXv2DQIJg9G/r1gyOOiDsiEUkStRSkftdfD0ceCWvWwN/+BuPGQfv2cUclIkkSxx7N+5rZ9ISv1WZ2pZntZmYTzWx+9KqprXFxh6qqcHzggTBkSGglnHFGvHGJSNKlPCm4+zx37+XuvYAiYC3wLHAdMMnduwOToveSah9/DGeeCXfdFd6fc044rmO9IhHJPnF3H50AfODuC4G+wJiofAxwVmxR5aKNG+Huu6FHD3j1Vdh++7gjEpEYxD3Q3B94PDpu5+5LAdx9qZntXtcNZjYIGATQqVOnlASZ9WbOhAsugLIyOP10uP9+6NIl7qhEJAaxtRTMrCVwJvDUltzn7qPcvdjdiwsKCpITXK75/HNYsgSefDIMJishiOSsOFsKpwHT3H1Z9H6ZmbWPWgntgeUxxpb9/v53mDYtzEY+4oiwZpFmJIvkvDjHFM5lU9cRwPPAwOh4IDA+5RHlgv/+F37wg/Ak0RNPwFdfhXIlBBEhpqRgZjsAJwHPJBSPAE4ys/nRuRFxxJa1qqpg1CjYf/+wN/Itt4QxBCUDEUkQS/eRu68F2tQqqyA8jSTJsHgxXHFF2BbzoYdgn33ijkhE0lDcj6RKMn31Ffz5z2EyWqdOMGVKeNxUCUFENkNJIVu9/jocfHDY72DKlFB24IHa+EZE6qWkkG0qKsKcg+OOgw0b4OWXwy5oIiKNEPfkNWlO7iEZzJkD114LN9wAO+wQd1QikkGUFDJc2cJVzP7XDHr26UlRtwK4805o1y5sfiMisoXUfZTByhYs59ULrqLfj09l4sW/omzhqrBfshKCiDSRkkKmmjKFzqccwzWvPso/ux7CC/scRWl5RdxRiUiGU/dRJrrzTrjmGnZutweXff/XvLR3CS3y8ygpbNPwvSIi9VBSyCSVlZCfD4cfDoMH03L4cM5ftZH9yysoKWxDUWftSyQiW0dJIRN88glcfjl06AD33RcWsIv2SS7aGSUDEWk2GlNIZxs3hr0N9t8fXnoJOneOOyIRyXJqKaSrefNg4EB4553wRNEDD0C3bnFHJSJZTkkhXbVoAUuXhrWLfvhDLU8hIimhpJBOJkyA8ePDuEFhIXzwQRhYFhFJEY0ppIPly+FHP4JTToFXXgnrF4ESgoikXFyb7OxqZk+b2ftmNtfMDjez3cxsopnNj16z/5Ead3j00TCQPG4cXH89zJgBbdvGHZmI5Ki4Wgr3AP9w9/2Ag4G5wHXAJHfvDkyK3me31ath6NCQFKZPD7uhaSc0EYlRypOCme0MfAt4GMDd17v7p0BfYEx02RjgrFTHlhLr14cniSorYZdd4K23YPJk6NEj7shERGJpKRQCK4BHzexdMxttZq2Adu6+FCB63T2G2JLrzTehVy+49FL4+99DWbdukKehHRFJD3H8NsoHegMPuPshwBdsQVeRmQ0ys6lmNnXFihXJirF5ffop/OxncPTRsHZtSAhnnhl3VCIi3xBHUlgMLHb3d6L3TxOSxDIzaw8QvS6v62Z3H+Xuxe5eXFBQkJKAt1q/fjB6NFx1FcyeDaedFndEIiJ1Svkzj+7+XzNbZGb7uvs84ARgTvQ1EBgRvY5PdWzNauFCaN0adt4ZRowIk8969447KhGResXVmX05MNbMZgK9gFsJyeAkM5sPnBS9zzyVlXDXXWHg+IYbQllRkRKCiGSEWGZHuft0oLiOUyekOpZmVVYGgwbBtGlwxhkwZEjcEYmIbBE99tJcHn4Y+vSBJUvCRLQXXtCqpiKScZQUtta6deH1uOPgkktg7twwsKwF7EQkAykpNNXSpXDOOXD22WG5isLCsJDdrrvGHZmISJMpKWypqip46KGwNMXzz0NJSSgTEckCWoZzSyxcCAMGhKUpjj02JId99ok7KhGRZqOksCV23hlWrQormw4cqHEDEck6Odt9VLZwFfe/toCyhavqv/C118LA8YYNYTLae+/BT36ihCAiWSknWwplC1cxYHQp6yuraJmfx9iLSijqXGv7hooKuPpqeOyxsGjdokVhMFmL14lIFsvJ33Cl5RWsr6yiymFDZRWl5RWbTrqHfZH32y+8Dh0aWgeFhfEFLCKSIjnZUigpbEPL/Dw2VFbRIj+PksI2m05WVoa1ivbeG0aNggMPjC9QEZEUy8mkUNS5NWMvKqG0vIKSwjYUddgR7r4bzj8/bHwzYQK0awfbbBN3qCIiKZWTSQFCYijq3BreeQfOHAQzZ0KrVvDTn0KHDnGHJyISi5wcUwDC/siXXQaHHw4rV8L48SEhiIjksNxNCpddBiNHwuWXw5w52glNRIQc7j7i5ptDYujTJ+5IRETSRu4mha5dw5eIiNSIJSmY2UfAGmAjUOnuxWa2G/Ak0AX4CDjH3RuYbiwiIs0pzjGF49y9l7tX78B2HTDJ3bsDk6L3IiKSQuk00NwXGBMdjwHOijEWEZGcFFdScGCCmZWZ2aCorJ27LwWIXnev60YzG2RmU81s6ooVK1IUrohIbohroPlId19iZrsDE83s/cbe6O6jgFEAxcXFnqwARURyUSwtBXdfEr0uB54F+gDLzKw9QPS6PI7YRERyWcqTgpm1MrOdqvcIxJ8AAAgrSURBVI+Bk4FZwPPAwOiygcD4VMcmIpLr4ug+agc8a2GTmnzgL+7+DzP7NzDOzC4EPgb6xRCbiEhOM/fM7ZY3sxXAwq34iLbA/5opnDhlSz1AdUlH2VIPUF2qdXb3grpOZHRS2FpmNjVhnkTGypZ6gOqSjrKlHqC6NEY6zVMQEZGYKSmIiEiNXE8Ko+IOoJlkSz1AdUlH2VIPUF0alNNjCiIi8nW53lIQEZEESgoiIlIja5OCmXU0s9fMbK6ZzTazK6Ly35jZTDObbmYTzKxDwj1DzWyBmc0zs1Pii36TzdUj4fzVZuZm1jahLO3qAfX+TG4ys0+in8l0Mzs94Z6Mqkt07vIo3tlmdltCeUbVxcyeTPiZfGRm0xPuSbu61FOPXmZWGtVjqpn1Sbgn7eoB9dblYDN728zeM7MXzGznhHuapy7unpVfQHugd3S8E/AfoAewc8I1PwcejI57ADOAbYGuwAfANulaj+h9R+BlwgS+tulcjwZ+JjcBV9dxfSbW5TjgFWDb6NzumVqXWtfcCdyQznWp52cyATgtKj8deD2d69FAXf4NHBOVXwD8prnrkrUtBXdf6u7TouM1wFxgT3dfnXBZK8Iy3hD2c3jC3de5+4fAAsJCfbHaXD2i078DfsmmOkCa1gMarEtdMrEulwAj3H1ddK56YcdMrAsAFtakOQd4PCpKy7rUUw8Hqv9HvQuwJDpOy3pAvXXZF5gcXTYRODs6bra6ZG1SSGRmXYBDgHei98PNbBEwALghumxPYFHCbYup/xdWyiXWw8zOBD5x9xm1Lkv7esA3fybAZVG33iNm1joqy8S67AMcbWbvmNkbZnZodFkm1qXa0cAyd58fvU/7utSqx5XA7dG/+TuAodFlaV8P+EZdZgFnRqf6EXoLoBnrkvVJwcx2BP4KXFndSnD3Ye7eERgLXFZ9aR23p83zuon1ACqBYWxKaF+7tI6ytKkH1PkzeQDoBvQClhK6KiAz65IPtAZKgGsIizwamVmXaueyqZUAaV6XOupxCTAk+jc/BHi4+tI6bk+bekCddbkAGGxmZYRupfXVl9Zxe5PqktVJwcxaEP5Ax7r7M3Vc8hc2Nb8WsynrAuzFpmZmrOqoRzdCv+EMM/uIEOs0M9uDNK4H1P0zcfdl7r7R3auAP7Cp2ZtxdSHE/IwHU4AqwsJlmVgXzCwf+B7wZMLlaVuXzdRjIFB9/BQZ/PfL3d9395PdvYiQqD+ILm++usQ9oJKsL0Lm/CNwd63y7gnHlwNPR8c9+fpATTlpMOi0uXrUuuYjNg00p2U9GviZtE84HkLoG83UulwM3BId70No0lsm1iU6dyrwRq2ytKxLPT+TucCx0fEJQFk616OBulQ/uJAXnb+guesSe+WT+Id6FKH5NBOYHn2dTsi8s6LyFwiDz9X3DCNk3nlETyvE/bW5etS6piYppGs9GviZ/Al4Lyp/vlaSyLS6tAT+HP0dmwYcn6l1ic49Blxcxz1pV5d6fiZHAWXRL813gKJ0rkcDdbmC8CTSf4ARRKtSNGddtMyFiIjUyOoxBRER2TJKCiIiUkNJQUREaigpiIhIDSUFERGpoaQgGcnM2iSs4PnfWqustow7vrqY2QXRBMNkfX4rM3vdzPLMbO9aq5pebGb/NrNdzOxuM/tWsuKQzJYfdwAiTeHuFYRlMTCzm4DP3f2OWIMKsWzj7hs3c/oCwtyF/27B5+W7e2UjL78IeMrdq8LKGjWfcT5hUt3x7v6Zmd0L3MemhdVEaqilIFnHzAaa2ZSo1TAy+p9zvpl9ama3m9k0M3vZzA6LFq0rt2gPBzO7yMyejc7PM7NfN/Jzf2tmU4A+ZnZz9L/yWWb2oAU/ICSx6j0KWprZYjPbNfrsEjN7JTr+rZk9ZGYTgUej73FX9L1nmtlFm6n6AGB8rT+LHwJXASe7+0oAd/8AaG9mBc36By9ZQUlBsoqZHQB8FzjC3XsRWsP9o9O7ABPcvTdhIbGbCMse9ANuSfiYPtE9vYEfWtikpaHPnebufdz9beAedz8UODA6d6q7P0mYlfoDd+/l7uup3yHAd9z9PGAQsNzd+wCHEhZE61Sr3tsBe7n74oTiQuAuQkJYzte9CxzRQAySg9R9JNnmRMIvzqlRF8r2bFpS+Et3nxgdvwd85u6VZvYe0CXhM15291UAZvYcYcmB/Ho+dz3wbML9J5jZNcB2hMXwyoCXtrAe4939q+j4ZGB/M0tMQt2BjxOu3x1YWeszlgFrCIs+3lvr3HKgAyK1KClItjHgEXe//muFYbXPxP+dVwHrEo4T/y3UXvvFG/jcLz1aL8bMdiD01/d290/M7LeE5FCXSja11mtf80WtOl3q7pM28zkAX27mM04D3jSz5VFrpdp20T0iX6PuI8k2rwDnWLRndfSUUqcG7qntZDPbNfoF3xd4aws+d3tCkvmfme3EpqXZIfyvfaeE9x8BRdFx4nW1vQxcGiUgzGxfM9s+8QJ3XwFsV/vJK3dfRljp9HYzOzHh1D6ERftEvkYtBckq7v6emd0MvGJmecAGwpM3W7K2/JuEvTa6AX9y9+kAjflcd68wszGEX7gL+foOZo8Co83sS8K4xU3AH8zsv8CUeuJ5COgETI+6rpYTklVtkwjjBK/XiukDMzsLeMHM+hK6zroQxhVEvkarpIokiJ7sOcDdr4w7li1lYevPS939/Aau6wf0cPebUxOZZBJ1H4lkCXf/N2H8oKF/1wb8LgUhSQZSS0FERGqopSAiIjWUFEREpIaSgoiI1FBSEBGRGkoKIiJS4/8DzXxERiArf3cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "p = np.polyfit(T,B,1)\n",
    "plot(T,B,'.', T,polyval(p,T),'r--')\n",
    "xlabel('Temperature (K)')\n",
    "ylabel('Radiance')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11a549f28>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfXRddZ3v8fc3Jz1pQcASmlhM2jQFy5OQkhADdwqj2IpzL0W8PPdeK9ByEWHQYbyji3sdGIclcx0YXchMNQWp3MqDgNBxSSlUkMwsYpu0RygypZimtIBJCUV6eTpJzvf+sXdOk5An2nOyz8PntVbWPmc/5ftLoJ/8fnuf3zZ3R0REBKAk6gJERCR3KBRERCRNoSAiImkKBRERSVMoiIhIWmnUBRyII444wmtqaqIuQ0Qkr7S3t7/u7jNG2pbXoVBTU0NbW1vUZYiI5BUz2zHaNg0fiYhImkJBRETSFAoiIpKW19cURtLb28uuXbt47733oi4l702dOpWqqiqmTJkSdSkiMkkKLhR27drFIYccQk1NDWYWdTl5y93p6elh165dzJkzJ+pyRGSSFNzw0XvvvUd5ebkC4QCZGeXl5epxiRSZggsFQIGQIfo5iuSm9h17uP3Jl2jfsSfj5y644SMRkULWvmMPS1a2kuxLES8tYfWyJupnT8/Y+QuypxC1N998k/POO49jjjmGY489lmeeeQaACy+8kLq6Ourq6qipqaGurg6Azs5Opk2blt525ZVXZr3Gyy67jIqKCk444YSsfy8RyZzWjh6SfSlSDr19KVo7ejJ6fvUUsuDaa6/lrLPO4oEHHiCZTPLOO+8AcN9996X3ue666zjssMPS7+fOnUsikZi0Gr/85S9z9dVX86UvfWnSvqeIHLim2nLipSX09qWYUlpCU215Rs+vUCDojrV29NBUW37A3bC33nqLp59+mrvuuguAeDxOPB4fso+7c//99/PrX//6Q51748aNXHvttbz99tuUlZWxfv16HnzwQR5++GH6+/vZsmUL1113HclkkrvvvpuysjJ+9atfcfjhh3/gXKeffjqdnZ3720wRiUj97OmsXtaUsX+zhiv64aOB8blb1m1lycrWA75w09HRwYwZM7j00kuZP38+y5Yt4+233x6yT0tLC5WVlRx99NHpddu3b2f+/PmcccYZtLS0fOC8yWSSCy+8kB/84Af87ne/44knnmDatGkAbNmyhZ/97Gds2LCB66+/noMOOojNmzdz6qmn8tOf/vSA2iMiuad+9nS++umjMh4IoFDI+PhcX18fmzZt4itf+QqbN2/m4IMP5uabbx6yzz333MPFF1+cfj9z5kxefvllNm/ezK233soll1zCW2+9NeSYrVu3MnPmTE455RQADj30UEpLg47epz/9aQ455BBmzJjBYYcdxtlnnw3AJz/5SfUGRArRzg3QckuwzLCiD4WB8bmYkZHxuaqqKqqqqvjUpz4FwHnnncemTZvS2/v6+njooYe48MIL0+vKysooLw++b319PXPnzuXFF18ccl53H/UW0bKysvTrkpKS9PuSkhL6+vrYuXNn+iL2ihUrDqh9IhKxnRtg1WL49U3BMsPBUPTXFDI9Pvexj32M6upqtm7dyrx581i/fj3HHXdcevsTTzzBMcccQ1VVVXrd7t27Ofzww4nFYnR0dLBt2zZqa2uHnPeYY47h1VdfZePGjZxyyins3bs3PXw0nurq6km9iC0iWdTZAv1J8P5g2dkC1Y0ZO33RhwIEwZDJsbnbbruNJUuWkEwmqa2t5Sc/+Ul627333jtk6Ajg6aef5tvf/jalpaXEYjFWrFjxgYvD8Xic++67j2uuuYZ3332XadOm8cQTT+x3jRdffDFPPfUUr7/+OlVVVdx4441cfvnl+30+EZkkNQsgFg8CIRYP3meQuXtGTziZGhoafPhDdl544QWOPfbYiCoqPPp5iuSgnRuCHkLNgv3qJZhZu7s3jLRNPQURkXxT3ZjRIaPBiv5Cs4iI7KNQEBHJM4nuBCufW0miO/M3kGj4SEQkjyS6Eyxft5xkf5J4LE7zombqKuoydn71FERE8khbVxvJ/iQpUvSmemnrahv/oA9BoSAikkcaKhuIx+LELMaUkik0VI54E9F+Uyhk2NatW9OfHq6rq+PQQw/l+9//PgBvvPEGCxcu5Oijj2bhwoXs2bNvnqXvfve7HHXUUcybN4/HHnss63Vq6myR/FRXUUfzomaunn91xoeOQKGQcfPmzSORSJBIJGhvb+eggw7i3HPPBeDmm2/mzDPPZNu2bZx55pnpOZF+//vfc++99/L888+zdu1arrrqKvr7+7Na55e//GXWrl2b1e8hItlRV1HHsk8uy3gggEIhkKXJpdavX8/cuXOZPXs2AI888ghLly4FYOnSpTz88MPp9RdddBFlZWXMmTOHo446ig0bPljLxo0bOe200zjppJNobGxk79693HXXXXzhC1/g7LPPZs6cOfzwhz/k1ltvZf78+TQ1NfHGG2+MWNvpp58+4pTaIlLcFApZnFxq+JQWXV1dzJw5EwhmRu3u7gbglVdeobq6Or1fVVUVr7zyypBzaepsEZkMCoWRJpfKgGQyyZo1azj//PPH3XekqUaGz4iqqbNFZDIoFAYml7JYRieXevTRRzn55JOprKxMr6usrOS1114D4LXXXqOiogIIegY7d+5M77dr1y6OPPLIIefT1NkiMhkUCtWNsHQNfOb6YJmh+USGP0gHYPHixaxatQqAVatWcc4556TX33vvvbz//vts376dbdu20dg4tI7BU2cD7N27l76+vgnVMjB1diKR4MorrzzQpolIAdMnmiHjk0u98847PP744/zoRz8asv6b3/wmF1xwAXfccQezZs3i5z//OQDHH388F1xwAccddxylpaXcfvvtxGKxIcdq6mwRmQyaOlvGpJ+nSOEZa+psDR+JiEiaQkFERNIUCiIikqZQEBGRtKyFgpndaWbdZrZl0Lrvmdl/mNmzZvYLM/vooG3fMrOXzGyrmX0uW3WJiMjostlTuAs4a9i6x4ET3P1E4EXgWwBmdhxwEXB8eMw/m1kMERGZVFkLBXd/Gnhj2Lp17j7wiatWoCp8fQ5wr7u/7+7bgZeA7DyVOsvGmjr7G9/4Bscccwwnnngi5557Lm+++SYAnZ2dTJs2LX3MZHzATFNni8hIorymcBnwaPj648DOQdt2hevyzlhTZy9cuJAtW7bw7LPP8olPfILvfve76ePmzp2bPm4ypqLQ1NkiMpJIPtFsZtcDfcDqgVUj7Dbip+rM7ArgCoBZs2ZlpJ5Ed4K2rjYaKhsyOj/58KmzFy1alN7W1NTEAw888KHOt3HjRq699lrefvttysrKWL9+PQ8++CAPP/ww/f39bNmyheuuu45kMsndd99NWVkZv/rVr0acIvv000/XZHki8gGT3lMws6XAfwGW+L6PU+8CqgftVgW8OtLx7v5jd29w94YZM2YccD0DD8G+bdNtLF+3nER34oDPOWD41NmD3XnnnXz+859Pv9++fTvz58/njDPOoKXlgzO1aupsEZkMkxoKZnYW8DfAYnd/Z9CmNcBFZlZmZnOAo4HMPvFmFNl6CPZYU2ffdNNNlJaWsmTJEiB4tsLLL7/M5s2bufXWW7nkkkt46623hhyjqbNFZDJk85bUe4BngHlmtsvMLgd+CBwCPG5mCTNbAeDuzwP3A78H1gJfdffsPo8ylK2HYI80dTYEs6P+8pe/ZPXq1empsMvKyigvLwegvr6euXPn8uKLLw45TlNnixyY9h17uP3Jl2jfsWf8nYtY1q4puPtI4yZ3jLH/TcBN2apnNAMPwc70NYWRps5eu3Yt//AP/8BvfvMbDjrooPT63bt3c/jhhxOLxejo6GDbtm3U1tYOOXbw1NmnnHIKe/fuTQ8fjWdg6myRYtW+Yw9LVraS7EsRLy1h9bIm6mdPj7qsnKRPNJP5h2APTJ39xS9+ccj6q6++mr1797Jw4cIht54+/fTTnHjiiZx00kmcd955rFix4gMXhwdPnX3SSSexcOFC3nvvvf2u8eKLL+bUU09l69atVFVVcccdo+a1SN5r7egh2Zci5dDbl6K1oyfqknKWps6WMennKYVgoKfQ25diinoKY06drYfsiEjBq589ndXLmmjt6KGptryoA2E8CgURKQr1s6crDCagIK8p5POQWC7Rz1EKys4N0HJLsJRRFVxPYerUqfT09FBeXj7qLZwyPnenp6eHqVOnRl2KyIHbuQFWLYb+JMTisHRNRp/LXkgKLhSqqqrYtWsXu3fvjrqUvDd16lSqqqrG31Ek13W2BIHg/cGys0WhMIqCC4UpU6YwZ86cqMsQkVxSsyDoIQz0FGoWRF1Rziq4UBAR+YDqxmDIqLMlCAT1EkalUBCR4lDdqDCYgIK8+0hERPaPQkFEikKiO8HK51ZmdHr8QqThIxEpeAPPTUn2J4nH4jQvas7oA7UKiXoKIlLwsvXclEKkUBCRgpet56YUIg0fiUjBy9ZzUwqRQkFEikJdRZ3CYAI0fCQiImkKBRERSVMoiIhImkJBRETSFAoiIpKmUBARkTSFgoiIpCkUREQkTaEgIiJpCgUREUlTKIiISJpCQURE0hQKIiKSplAQEZE0hYKIiKQpFERkVO079nD7ky/RvmNP1KXIJNFDdkRkRO079rBkZSvJvhTx0hJWL2uifvb0qMuSLFNPQURG1NrRQ7IvRcqhty9Fa0dP1CXJJFAoiMiImmrLiZeWEDOYUlpCU2151CXJJNDwkYiMqH72dFYva6K1o4em2nINHRUJhYKIjKp+9nSFQZHJ2vCRmd1pZt1mtmXQusPN7HEz2xYupw/a9i0ze8nMtprZ57JVl4iIjC6b1xTuAs4atu6bwHp3PxpYH77HzI4DLgKOD4/5ZzOLZbE2EZmInRug5ZZgKUUha8NH7v60mdUMW30O8Ofh61XAU8DfhOvvdff3ge1m9hLQCDyTrfpEZBw7N8CqxdCfhFgclq6B6saoq5Ism+y7jyrd/TWAcFkRrv84sHPQfrvCdSISlc6WIBC8P1h2tkRdkUyCCYWCmVWa2R1m9mj4/jgzuzyDddgI63yUWq4wszYza9u9e3cGSxCRIWoWBD0EiwXLmgVRVySTYKI9hbuAx4Ajw/cvAl/bj+/XZWYzAcJld7h+F1A9aL8q4NWRTuDuP3b3BndvmDFjxn6UICITUt0YDBl95noNHRWRiYbCEe5+P5ACcPc+oH8/vt8aYGn4einwyKD1F5lZmZnNAY4GdGVLJGrVjbDgOgVCEZnohea3zayccEjHzJqAP411gJndQ3BR+Qgz2wX8LXAzcH849PQycD6Auz9vZvcDvwf6gK+6+/6EjoiIHICJhsJfEfw1P9fM/h2YAZw31gHufvEom84cZf+bgJsmWI+ITIJEd4K2rjYaKhuoq6iLuhyZBBMKBXffZGZnAPMILgpvdfferFYmIpFKdCdYvm45yf4k8Vic5kXNCoYiMNG7j74KfMTdn3f3LcBHzOyq7JYmIlFq62oj2Z8kRYreVC9tXW1RlySTYKIXmpe7+5sDb9x9D7A8OyWJSC5oqGwgHosTsxhTSqbQUNkQdUkyCSZ6TaHEzMzdBy40x4B49soSkajVVdTRvKhZ1xSKzERD4TGCu4ZWENyBdCWwNmtViUhOqKuoUxgUmYmGwt8A/wP4CsGF5nXAymwVJSIi0Zjo3Ucp4F/CLxERKVATCgUz+0/ADcDs8BgD3N1rs1eaiIhMtokOH90BfB1oZ/+mtxARkTww0VD4k7s/mtVKREQkchMNhSfN7HvAQ8D7AyvdfVNWqhIRkUhMNBQ+FS4Hf3rFgc9kthwREYnSRO8++nS2CxERkehN+BnNZvafgeOBqQPr3P3vslGUiIhEY6IT4q0ALgSuIbgd9XyC21NFRKSATHRCvNPc/UvAHne/ETiVoY/PFBGRAjDRUHg3XL5jZkcCvcCc7JQkIiJRmeg1hV+a2UeB7wGbCO480txHIqNo37GH1o4emmrLqZ89PepyRCZsoncffSd8+aCZ/RKY6u5jPqNZpFi179jDkpWtJPtSxEtLWL2sScEgeWPMUDCzz7j7r83siyNsw90fyl5pIvmptaOHZF+KlENvX4rWjh6FguSN8XoKZwC/Bs4eYZsTfMJZRAZpqi0nXlpCb1+KKaUlNNWWR12SyIRZ+DC1vNTQ0OBtbXpurOQeXVOQXGZm7e4+4vNVxxs++quxtrv7rQdSmEihqp89XWEgeWm84aNDwuU84BRgTfj+bODpbBUlIiLRGDMUwg+qYWbrgJPdfW/4/gbg51mvTiRf7dwAnS1QswCqG6OuRmTCJvo5hVlActD7JFCT8WpECsHODbBqMfQnIRaHpWsUDJI3JhoKdwMbzOwXBHcdnQv8NGtVieSzzpYgELw/WHa2KBQkb0z0w2s3mdla4M/CVZe6++bslSWSx2oWBD2EgZ5CzYKoKxKZsAlPne3u7Wa2k3DqbDOb5e4vZ60ykXxV3RgMGemaguShCYWCmS0GbgGOBLoJrjH8B8HzFURkuOpGhYHkpYnOkvodoAl40d3nAJ8F/j1rVYmISCQmGgq97t4DlJhZibs/CdRlsS6RvJboTrDyuZUkuhNRlyLyoUz0msKbZvYRgg+srTazbqAve2WJ5K9Ed4Ll65aT7E8Sj8VpXtRMXYX+hpL8MNGewjnAO8DXgbXAHxh5kjyRotfW1UayP0mKFL2pXtq6ND+X5I+J3pL6dvgyBawysxhwEbA6W4WJ5KuGygbisTi9qV6mlEyhoXLEecdEctJ4E+IdCnwV+DjBvEePh++/ASRQKIh8QF1FHc2LmmnraqOhskFDR5JXxusp3A3sAZ4BlhGEQRw4x911BU1kFHUVdQoDyUvjhUKtu38SwMxWAq8DswYmxhMRkcIy3oXm3oEX7t4PbM9EIJjZ183seTPbYmb3mNlUMzvczB43s23hUpPRi4hMsvFC4SQzeyv82gucOPDazN7an29oZh8H/hJocPcTgIGL1t8E1rv70cD68L2IiEyiMUPB3WPufmj4dYi7lw56fegBfN9SYJqZlQIHAa8S3Pa6Kty+CvjCAZxfRET2w0Q/p5Ax7v4K8I/Ay8BrwJ/cfR1Q6e6vhfu8BlRMdm0iIsVu0kMhvFZwDjCHYIK9g83sv32I468wszYza9u9e3e2yhQRKUqTHgoEk+ltd/fd7t4LPAScBnSZ2UyAcNk90sHu/mN3b3D3hhkzZkxa0SIixSCKUHgZaDKzg8zMgDOBFwg+HLc03Gcp8EgEtYmIFLUJP2QnU9z9t2b2ALCJYFK9zcCPgY8A95vZ5QTBcf5k1yYiUuwmPRQA3P1vgb8dtvp9gl6DiIhEJIrhI5ERte/Yw+1PvkT7jj1RlyJStCLpKYgM175jD0tWtpLsSxEvLWH1sibqZ+tD7SKTTT0FyQmtHT0k+1KkHHr7UrR29ERdkkhRUihITmiqLSdeWkLMYEppCU215VGXJFKUNHwkOaF+9nRWL2uitaOHptpyDR2JREShIDmjfvZ0hYFIxDR8JCIiaQoFyR07N0DLLcFSRCKh4SPJDTs3wKrF0J+EWByWroHqxqirEik66ilIbuhsCQLB+4NlZ0vUFYkUJYWC5IaaBUEPwWLBsmZB1BWJFCUNH0luqG4Mhow6W4JA0NCRSCQUCpI7qhsVBiIR0/CRiIikKRQkZyS6E6x8biWJ7kTUpYgULQ0fSU5IdCdYvm45yf4k8Vic5kXN1FXURV2WSNFRT0FyQltXG8n+JClS9KZ6aetqi7okkaKkUJCc0FDZQDwWJ2YxppRMoaGyIeqSRIqSho8kJ9RV1NG8qJm2rjYaKhs0dCQSEYWC5Iy6ijqFgUjENHwkIiJpCgUREUlTKIiISJpCQURE0hQKIiKSplAQEZE0hYKIiKQpFPJc+4493P7kS7Tv2BN1KSJSAPThtTzWvmMPS1a2kuxLES8tYfWyJupnT4+6LBHJY+op5LHWjh6SfSlSDr19KVo7eqIuSUTynEIhjzXVlhMvLSFmMKW0hKba8qhLEpE8p+GjPFY/ezqrlzXR2tFDU225ho5E5IApFPJc/ezpCgMRyRgNH4mISJpCId/t3AAttwRLEZEDpOGjfLZzA6xaDP1JiMVh6Rqoboy6KhHJY+op5LPOliAQvD9YdrZEXZGI5DmFQj6rWRD0ECwWLGsWRF2RiOS5SIaPzOyjwErgBMCBy4CtwH1ADdAJXODumrthLNWNwZBRZ0sQCBo6EpEDFFVP4QfAWnc/BjgJeAH4JrDe3Y8G1ofvZTzVjbDgOgWCiGTEpIeCmR0KnA7cAeDuSXd/EzgHWBXutgr4wmTXJiJS7KLoKdQCu4GfmNlmM1tpZgcDle7+GkC4rIigtryT6E6w8rmVJLoTUZciIgUgimsKpcDJwDXu/lsz+wEfYqjIzK4ArgCYNWtWdirME4nuBMvXLSfZnyQei9O8qJm6irqoyxKRPBZFT2EXsMvdfxu+f4AgJLrMbCZAuOwe6WB3/7G7N7h7w4wZMyal4FzV1tVGsj9JihS9qV7autqiLklE8tykh4K7/xHYaWbzwlVnAr8H1gBLw3VLgUcmu7Z801DZQDwWJ2YxppRMoaGyIeqSRCTPRfWJ5muA1WYWBzqASwkC6n4zuxx4GTg/otryRl1FHc2LmmnraqOhskFDRyJywCIJBXdPACP9WXvmZNeS7+oq6hQGIpIx+kSziIikKRRERCRNoSAiImkKBRERSVMoiIhIWtGGQvuOPdz+5Eu079BErCIiA4ryyWvtO/awZGUryb4U8dISVi9ron729KjLEhGJXFH2FFo7ekj2pUg59PalaO3oibokEZGcUJSh0FRbTry0hJjBlNISmmrLoy5JRCQnFOXwUf3s6axe1kRrRw9NteUaOhIRCRVlKEAQDAoDEZGhinL4SERERla8obBzA7TcEixFRAQo1uGjnRtg1WLoT0IsDkvX6MH3IiIUa0+hsyUIBO8Plp0tUVckIpITijMUahYEPQSLBcuaBVFXJCKSE4pz+Ki6MRgy6mwJAkFDRyIiQLGGAgRBoDAQERmiOIePgER3gpXPrSTRnYi6FBGRnFGUPYVEd4Ll65aT7E8Sj8VpXtSs5xyLiFCkPYW2rjaS/UlSpOhN9dLW1RZ1SSIiOaEoQ6GhsoF4LE7MYkwpmUJDZUPUJYmI5ISiHD6qq6ijeVEzbV1tNFQ2aOhIRCRUlKEAQTAoDEREhirK4SMRERmZQkFERNIUCiIikqZQEBGRNIWCiIikKRRERCTN3D3qGvabme0GdhzAKY4AXs9QOVEqlHaA2pKLCqUdoLYMmO3uM0bakNehcKDMrM3d8/7jzIXSDlBbclGhtAPUlonQ8JGIiKQpFEREJK3YQ+HHUReQIYXSDlBbclGhtAPUlnEV9TUFEREZqth7CiIiMohCQURE0go2FMys2syeNLMXzOx5M7s2XP8dM3vWzBJmts7Mjhx0zLfM7CUz22pmn4uu+n1Ga8eg7X9tZm5mRwxal3PtgDF/JzeY2Svh7yRhZn8x6Ji8aku47Zqw3ufN7P8MWp9XbTGz+wb9TjrNLDHomJxryxjtqDOz1rAdbWbWOOiYnGsHjNmWk8zsGTN7zsz+1cwOHXRMZtri7gX5BcwETg5fHwK8CBwHHDpon78EVoSvjwN+B5QBc4A/ALFcbUf4vhp4jOADfEfkcjvG+Z3cAPz1CPvnY1s+DTwBlIXbKvK1LcP2uQX4di63ZYzfyTrg8+H6vwCeyuV2jNOWjcAZ4frLgO9kui0F21Nw99fcfVP4ei/wAvBxd39r0G4HAwNX2s8B7nX39919O/AS0EjERmtHuPmfgP/JvjZAjrYDxm3LSPKxLV8Bbnb398Nt3eEh+dgWAMzMgAuAe8JVOdmWMdrhwMBf1IcBr4avc7IdMGZb5gFPh7s9DvzX8HXG2lKwoTCYmdUA84Hfhu9vMrOdwBLg2+FuHwd2DjpsF2P/gzXpBrfDzBYDr7j774btlvPtgA/+ToCrw2G9O81serguH9vyCWCBmf3WzH5jZqeEu+VjWwYsALrcfVv4PufbMqwdXwO+F/4//4/At8Ldcr4d8IG2bAEWh5vOJxgtgAy2peBDwcw+AjwIfG2gl+Du17t7NbAauHpg1xEOz5n7dQe3A+gDrmdfoA3ZdYR1OdMOGPF38i/AXKAOeI1gqALysy2lwHSgCfgGcH/4l3Y+tmXAxezrJUCOt2WEdnwF+Hr4//zXgTsGdh3h8JxpB4zYlsuAr5pZO8GwUnJg1xEO36+2FHQomNkUgh/oand/aIRdfsa+7tcu9qUuQBX7upmRGqEdcwnGDX9nZp0EtW4ys4+Rw+2AkX8n7t7l7v3ungKa2dftzbu2ENT8kAc2ACmCicvysS2YWSnwReC+QbvnbFtGacdSYOD1z8nj/77c/T/cfZG71xME9R/C3TPXlqgvqGTriyA5fwp8f9j6owe9vgZ4IHx9PEMv1HSQAxedRmvHsH062XehOSfbMc7vZOag118nGBvN17ZcCfxd+PoTBF16y8e2hNvOAn4zbF1OtmWM38kLwJ+Hr88E2nO5HeO0ZeDGhZJw+2WZbkvkjc/iD/XPCLpPzwKJ8OsvCJJ3S7j+XwkuPg8ccz1B8m4lvFsh6q/R2jFsn3Qo5Go7xvmd3A08F65fMywk8q0tceD/hv+NbQI+k69tCbfdBVw5wjE515Yxfid/BrSH/2j+FqjP5XaM05ZrCe5EehG4mXBWiky2RdNciIhIWkFfUxARkQ9HoSAiImkKBRERSVMoiIhImkJBRETSFAqSl8ysfNAMnn8cNstqPOr6RmJml4UfMMzW+Q82s6fMrMTMjho2q+mVZrbRzA4zs++b2enZqkPyW2nUBYjsD3fvIZgWAzO7Afh/7v6PkRYV1BJz9/5RNl9G8NmFP36I85W6e98Ed18G/NzdU8HMGulzXErwobrPuPufzOw24Ifsm1hNJE09BSk4ZrbUzDaEvYZ/Dv9yLjWzN83se2a2ycweM7NPhZPWdVj4DAczW2Zmvwi3bzWz/zXB8/69mW0AGs3sxvCv8i1mtnZu6REAAALFSURBVMICFxKE2MAzCuJmtsvMPhqeu8nMnghf/72Z/cjMHgd+En6PW8Pv/ayZLRul6UuAR4b9LC4BrgMWufsbAO7+B2Cmmc3I6A9eCoJCQQqKmZ0AnAuc5u51BL3hi8LNhwHr3P1kgonEbiCY9uB84O8GnaYxPOZk4BILHtIy3nk3uXujuz8D/MDdTwE+GW47y93vI/hU6oXuXufuScY2Hzjb3f87cAXQ7e6NwCkEE6LNGtbuqUCVu+8atLoWuJUgELoZajNw2jg1SBHS8JEUms8S/MPZFg6hTGPflMLvuvvj4evngD+5e5+ZPQfUDDrHY+6+B8DMHiaYcqB0jPMmgV8MOv5MM/sGMJVgMrx24NEP2Y5H3P298PUi4FgzGxxCRwMvD9q/Anhj2Dm6gL0Ekz7eNmxbN3AkIsMoFKTQGHCnu//vISuD2T4H/3WeAt4f9Hrw/wvD537xcc77rofzxZjZQQTj9Se7+ytm9vcE4TCSPvb11ofv8/awNl3l7utHOQ/Au6Oc4/PAv5lZd9hbGTA1PEZkCA0fSaF5ArjAwmdWh3cpzRrnmOEWmdlHw3/gzwH+/UOcdxpByLxuZoewb2p2CP5qP2TQ+06gPnw9eL/hHgOuCgMIM5tnZtMG7+Duu4Gpw++8cvcugplOv2dmnx206RMEk/aJDKGeghQUd3/OzG4EnjCzEqCX4M6bDzO3/L8RPGtjLnC3uycAJnJed+8xs1UE/+DuYOgTzH4CrDSzdwmuW9wANJvZH4ENY9TzI2AWkAiHrroJwmq49QTXCZ4aVtMfzOwLwL+a2TkEQ2c1BNcVRIbQLKkig4R39pzg7l+LupYPy4JHf17l7peOs9/5wHHufuPkVCb5RMNHIgXC3TcSXD8Y7/9rA/5pEkqSPKSegoiIpKmnICIiaQoFERFJUyiIiEiaQkFERNIUCiIikvb/AWyrvpyapaD0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure()\n",
    "plot(T, bt2rad(675., T), '.', T, bt2rad(700., T), '.', T, bt2rad(725., T), '.')\n",
    "xlabel('Temperature (K)')\n",
    "ylabel('Radiance')\n",
    "legend(('675 cm-1', '700 cm-1', '725 cm-1'), loc='best')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "p1 = np.polyfit(T, bt2rad(675., T), 1)\n",
    "p2 = np.polyfit(T, bt2rad(700., T), 1)\n",
    "p3 = np.polyfit(T, bt2rad(725., T), 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[array([   1.31488375, -250.3667531 ]),\n",
       " array([   1.31639086, -253.49107034]),\n",
       " array([   1.3123575 , -255.37736595])]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[p1, p2, p3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[280.02353471,   0.        ,   0.        ],\n",
       "       [  0.        , 279.95304433,   0.        ],\n",
       "       [  0.        ,   0.        , 280.62054301]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = 280.\n",
    "a = np.array([p1[0], p2[0], p3[0]])\n",
    "b = np.array([p1[1], p2[1], p3[1]])\n",
    "\n",
    "B = bt2rad(np.array([675., 700., 725.]), T) + randn(3,1)*0.5\n",
    "\n",
    "np.linalg.inv(a*np.identity(3)) * (B - b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
